/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/*********************************************************************************************************************
 * nebrList.c -- This module collects functions and variables related to the building of the neighbor list.          *
 *               A substantial part of the code has been adapted from 'The Art of Molecular Dynamics Simulation' by  *
 *               D.C. Rapaport (2004, 2nd ed, pp. 49-58).                                                            *  
 *********************************************************************************************************************/

#include "dataStructs.h"

// global variables
int *cellList; // pointer array: the first nSites elements are pointers linking the atoms belonging to the same cells, 
               // the remaining VProd(cells) elements, one per cell, point to the initial atom in the list for each cell.   
int *interactionType, nebrTabLen;
real *dispBetweenNebrListRefreshes; /* accumulated displacements of each site between neighbor-list reconstructions */
static int rebuildNebrListNow; /* Boolean: initialised to 0, it is set to 1 when nebrList refresh is due */ 

/* external data reference(s) */
extern const VecR *siteCoordinates; 
extern const VecR *siteVelocity; 
extern const Site *site;  /* array of Site, represents the state of the system */
extern const VecR region; /* edge lengths of the simulation region */ 
extern const real deltaT, rCutMax;
extern const int nSites, nebrTabMax, maxEdgeCells, stepCount;
extern VecI cells; // vector of integers indicating the number of cells along each dimension 
                   // (the product cells.x * cells.y * cells.z gives the total number of cells)
extern const real rNebrShell;
extern int *interactionType, *nebrTab, nebrTabLen;

void ErrExit(); // external function declaration

void UpdateCellSize() // [Rapaport, p. 161] 
{ 
  VSCopy( cells, 1. / rCutMax, region );
  cells.x = Min( cells.x, maxEdgeCells );
  cells.y = Min( cells.y, maxEdgeCells );
  cells.z = Min( cells.z, maxEdgeCells );
}

/*********************************************************************************************************************
 * DetermineInteractionType -- given the two types of a pair of sites, returns the type of nonbonded interaction 
 *********************************************************************************************************************/
static InteractionType DetermineInteractionType( SiteType typeSite1, SiteType typeSite2 ) 
{
  int typeSum = typeSite1 + typeSite2; /* used to identify the interaction type [Rapaport, p.219] */

  switch ( typeSum ) {
  case 0:  /* typeSum=0 implies  SSD (type 0) + SSD (type 0) pair */
    return SSD_SSD;  /* printf("## typeSite1 = %d, typeSite2 = %d (0 0)\n", typeSite1, typeSite2 ); */
    break;
  case 1: /* SSD water (type 0) + CHOL (type 1) */
    if (typeSite1 == WATER_TYPE) return DIP_CHG; /* water-choline */
    else return CHG_DIP;  /* choline-water printf("## typeSite1 = %d, typeSite2 = %d (1 0 or 0 1)\n", typeSite1, typeSite2 ); */
    break;
  case 2: /* two possibilities: typeSum = 0 (SSD) + 2 (PHOS) or typeSum = 1 (CHOL) + 1 (CHOL) */
    if ( (typeSite1 == CHOLINE_TYPE) || (typeSite2 == CHOLINE_TYPE) ) { /* must be a choline-choline interaction */
      return CHG_CHG; /* printf("## typeSite1 = %d, typeSite2 = %d (1 1)\n", typeSite1, typeSite2 ); */
      break;
    } else { /* SSD (type 0) + PHOS (type 2) */
      if (typeSite1 == WATER_TYPE) return DIP_CHG; /* water-phosphate */
      else return CHG_DIP;  /* phosphate-water printf("## typeSite1 = %d, typeSite2 = %d (either 0 2 or 2 0)\n", typeSite1, typeSite2 ); */
      break;
    }
  case 3: /* two possibilities: either 0+3 or 1+2 */
    if ( (typeSite1 == WATER_TYPE) || (typeSite2 == WATER_TYPE) ) { 
      return DIP_DIP; /* if one is WATER (type 0) the other must be GLYCEROL (type 3) */ 
      break;
    } else { /* CHOL (type 1) + PHOS (type 2) */
      return CHG_CHG; /* printf("## typeSite1 = %d, typeSite2 = %d (either 1 2 or 2 1)\n", typeSite1, typeSite2 ); */
      break;
    }
  case 4: /* three possibilities: 0+4 || 1+3 || 2+2 */
    if ( (typeSite1 == WATER_TYPE) || (typeSite2 == WATER_TYPE) ) { /* WATER (type 0) + ESTER (type 4) */
      return DIP_DIP; /* printf("## typeSite1 = %d, typeSite2 = %d (either 0 4 or 4 0)\n", typeSite1, typeSite2 ); */
      break;
    } else if ( (typeSite1 == CHOLINE_TYPE) || (typeSite2 == CHOLINE_TYPE) ) { /* CHOL (type 1) + GLYCEROL (type 3) */
      if ( typeSite1 == CHOLINE_TYPE ) return CHG_DIP; 
      else return DIP_CHG; 
      break;
    } else { /* PHOS-PHOS (2+2) */
      return CHG_CHG; 	/* printf("## typeSite1 = %d, typeSite2 = %d (2 2)\n", typeSite1, typeSite2 ); */
      break;
    }
  case 5: // three possibilities 0+5 || 1+4 || 2+3
    if ( (typeSite1 == TAIL_TYPE) || (typeSite2 == TAIL_TYPE) ) { // WATER (type 0) + TAIL (type 5)
      return NEUTRAL; // printf("## typeSite1 = %d, typeSite2 = %d (either 0 5 or 5 0)\n", typeSite1, typeSite2 );
      break;
    } else { // CHOL (type 1) + ESTER (type 4) or PHOS (type2) + GLYCEROL (type 3)
      if ( (typeSite1 == CHOLINE_TYPE) || (typeSite1 == PHOSPHATE_TYPE) ) 
	return CHG_DIP; // printf("## typeSite1 = %d, typeSite2 = %d (either 3 2, 2 3 or 1 4, 4 1)\n", typeSite1, typeSite2 );
      else return DIP_CHG; 
      break;
    }
  case 6: // four possibilities: i) 1+5 (or 5+1) ii) 2+4 iii) 3+3 iv) 6+0
    if ( (typeSite1 == TAIL_TYPE) || (typeSite2 == TAIL_TYPE) ) { // CHOL (type 1) + TAIL (type 5)
      return NEUTRAL; // printf("## typeSite1 = %d, typeSite2 = %d (1 5 or 5 1)\n", typeSite1, typeSite2 );
      break;
    } else if ( (typeSite1 == PHOSPHATE_TYPE) || (typeSite2 == PHOSPHATE_TYPE) ) { // PHOS (type 2) + ESTER (type 4)
      if (typeSite1 == PHOSPHATE_TYPE) return CHG_DIP; // printf("## typeSite1 = %d, typeSite2 = %d (either 4 2 or 2 4)\n", typeSite1, typeSite2 );
      else return DIP_CHG; 
      break;
    } else if ( (typeSite1 == AMINE_TYPE) || (typeSite2 == AMINE_TYPE) ) { // AMINE (type 6) + SSD (type 0)
      if (typeSite1 == AMINE_TYPE) return CHG_DIP; // printf("## typeSite1 = %d, typeSite2 = %d (either 4 2 or 2 4)\n", typeSite1, typeSite2 );
      else return DIP_CHG; 
      break;
    } else { // GLYCEROL (type 3) + GLYCEROL (type 3)
      return DIP_DIP; // printf("## typeSite1 = %d, typeSite2 = %d (3 3)\n", typeSite1, typeSite2 );
      break;
    }
  case 7: // three possibilities: i) 2+5 (or 5+2) ii) 3+4 (or 4+3) iii) 6+1
    if ( (typeSite1 == PHOSPHATE_TYPE) || (typeSite2 == PHOSPHATE_TYPE) ) { // if one is PHOSPHATE (type 2) the other must be TAIL (type 5)
      return NEUTRAL; // printf("## typeSite1 = %d, typeSite2 = %d (2 5 or 5 2)\n", typeSite1, typeSite2 );
      break;
    } else if ( (typeSite1 == AMINE_TYPE) || (typeSite2 == AMINE_TYPE) ) { /* if either of the two is a AMINE (type 6), must be a choline-amine interaction */
      return CHG_CHG; /* printf("## typeSite1 = %d, typeSite2 = %d (1 1)\n", typeSite1, typeSite2 ); */
      break;
    } else { // GLYCEROL (type 3) + ESTER (type 4)
      return DIP_DIP; // printf("## typeSite1 = %d, typeSite2 = %d (3 4 or 4 3)\n", typeSite1, typeSite2 );
      break;
    }
  case 8: // three possibilities: i) 3+5 (or 5+3) ii) 4+4 iii) 6+2
    if ( (typeSite1 == ESTER_TYPE) || (typeSite2 == ESTER_TYPE) ) { // ESTER (type 4) + ESTER (type 4)
      return DIP_DIP; // printf("## typeSite1 = %d, typeSite2 = %d (4 4)\n", typeSite1, typeSite2 );
      break;
    } else if ( (typeSite1 == AMINE_TYPE) || (typeSite2 == AMINE_TYPE) ) { /* if either of the two is a AMINE (type 6), must be a phosphate-amine interaction */
      return CHG_CHG; /* printf("## typeSite1 = %d, typeSite2 = %d (1 1)\n", typeSite1, typeSite2 ); */
      break;
    } else { // GLYCEROL (type 3) + TAIL (type 5)
      return NEUTRAL; // printf("## typeSite1 = %d, typeSite2 = %d (3 5 or 5 3)\n", typeSite1, typeSite2 );
      break;
    }
  case 9: // two possibilities: i) 5 (TAIL) + 4 (ESTER) ii) 6 (AMINE) + 3 (GLYCEROL)
  case 10: // 5 (TAIL) + 5 (TAIL) or 6 (AMINE) + 4 (ESTER) 
    if ( (typeSite1 == AMINE_TYPE) || (typeSite2 == AMINE_TYPE) ) { // AMINE + GLYCEROL or AMINE + ESTER (charge-dipole pair in both cases)
      if (typeSite1 == AMINE_TYPE) return CHG_DIP; // printf("## typeSite1 = %d, typeSite2 = %d (either 4 2 or 2 4)\n", typeSite1, typeSite2 );
      else return DIP_CHG; 
      break;
    } else { // TAIL (type 5) + ESTER (4) or TAIL (type 5) + TAIL (5)
      return NEUTRAL; // printf("## typeSite1 = %d, typeSite2 = %d (3 5 or 5 3)\n", typeSite1, typeSite2 );
      break;
    }
  case 11: // AMINE - TAIL
    return NEUTRAL; // printf("## typeSite1 = %d, typeSite2 = %d (3 5 or 5 3)\n", typeSite1, typeSite2 );
    break;
  case 12: // AMINE - AMINE
    return CHG_CHG; 
    break;
  default: 
    printf("DetermineInteractionType() - typeSite1 = %d, typeSite2 = %d\n", typeSite1, typeSite2 );
    ErrExit( ERR_UNDETERMINED_INTERACTION ); 
  }  
  return -1; /* only for debugging purposes as execution should never get to this point */
}

/*********************************************************************************************************************
 * BuildNeighborList() -- Adapted from [Rapaport, p.57] 
 *********************************************************************************************************************/
static void BuildNeighborList() { 
  VecR dr, invWid, rs, shift;
  VecI cc, m1v, m2v, vOff[] = OFFSET_VALS;
  real rrNebr;
  int c, j1, j2, m1, m1x, m1y, m1z, m2, n, offset;

  rrNebr = Sqr( rCutMax + rNebrShell );
  VDiv( invWid, cells, region );  // invWid = 1 / cell-width
  for ( n = nSites; n < nSites + VProd( cells ); n++ ) cellList[ n ] = -1; // Initialization; VProd( cells ) is the total number of cells
  DO_SITE {
    VSAdd( rs, siteCoordinates[ n ], 0.5, region ); // rs := r[n] + edge/2, shifting coordinates to positive values 
    VMul( cc, rs, invWid ); // cc contains the (integer) "coordinates" identifying the cell containing current site
    c = VLinear( cc, cells ) + nSites; // the integer vector cc is transformed into a single integer, still identifying the cell containing current site
    cellList[ n ] = cellList[ c ]; // the first time a cell is considered cellList[n] = cellList[c] = (-1), hence n will be the last element of the cell-list;
    // possible extra elements are then added on to the most recently one added, ie, the list is constructed in a backward fashion
    cellList[ c ] = n;             // cellList[c] -> cellList[n] the first element of cell "c" is the one with the index "n";
    // this will change every time a new element is added to the list of the particular cell
  }
  nebrTabLen = 0;
  for ( m1z = 0; m1z < cells.z; m1z++ ) { // scan over cells
    for ( m1y = 0; m1y < cells.y; m1y++ ) {
      for ( m1x = 0; m1x < cells.x; m1x++ ) {
	VSet (m1v, m1x, m1y, m1z);
	m1 = VLinear (m1v, cells) + nSites;
	for (offset = 0; offset < N_OFFSET; offset++) { // scan over offsets
	  VAdd( m2v, m1v, vOff[ offset ] );
	  VZero( shift );
	  VCellWrapAll();
	  m2 = VLinear (m2v, cells) + nSites;
	  DO_CELL (j1, m1) { /* loop over all sites belonging to cell m1 */
	    DO_CELL (j2, m2) { /* loop over all sites belonging to cell m2 */
	      if ((m1 != m2 || j2 < j1) && 
		  /* exclude intralipid nearest neighbors (bonded pairs in the same lipid) [Rapaport, p.247]: */
		  /* proceed to consider pair interaction only if at least one of the following 3 conditions is true: */
		  ((!site[j1].inLipid) /* 1) j1 is not a lipid site */
		   || (site[j1].inLipid != site[j2].inLipid) /* 2) j1 and j2 belong to different lipids */
		   || ((site[j1].lipidUnit!=site[j2].skipNonbondedInteraction1)&& /* 3) site-specific exclusions */  
		       (site[j1].lipidUnit!=site[j2].skipNonbondedInteraction2)&& 
		       (site[j1].lipidUnit!=site[j2].skipNonbondedInteraction3)&&
		       (site[j2].lipidUnit!=site[j1].skipNonbondedInteraction1)&& 
		       (site[j2].lipidUnit!=site[j1].skipNonbondedInteraction2)&& 
		       (site[j2].lipidUnit!=site[j1].skipNonbondedInteraction3)))) { /* 'if' block: */
		VSub( dr, siteCoordinates[j1], siteCoordinates[j2] );
		VVSub( dr, shift );
		if ( VLenSq( dr ) < rrNebr ) {
		  if ( nebrTabLen >= nebrTabMax ) 
		    ErrExit( ERR_TOO_MANY_NEBRS );
		  nebrTab[ 2 * nebrTabLen ] = j1;
		  nebrTab[ 2 * nebrTabLen + 1 ] = j2;
		  interactionType[ nebrTabLen ] = DetermineInteractionType( site[j1].type, site[j2].type ); 
		  ++nebrTabLen;
		}
	      } /* end 'if' block: */
	    }
	  }
	}
      }
    }
  }
}

void InitNebrList() 
{
  int n;
  rebuildNebrListNow = 1;
  DO_SITE dispBetweenNebrListRefreshes[ n ] = 0; 
}


/************************************************************************************************************
 * FindLargest -- finds and returns the largest element in an array.
 ************************************************************************************************************/
real FindLargest(  int *indexLargestElement, real array[], int nElements ) 
{
  int i; // counter
  real largest = 0;
  
  for ( i = 0; i < nElements; ++i ) {
    if ( array[ i ] > largest ) {
      largest = array[ i ];
      *indexLargestElement = i;
    }
  }
  return largest;
}

/************************************************************************************************************
 * FindSecondLargest -- finds and returns the "second largest" element in an array.
 ************************************************************************************************************/
real FindSecondLargest( int indexLargestElement, real array[], int nElements )
{
  int i; // counter
  real secondLargest = 0;

  for ( i = 0; i < nElements; ++i ) {
    if ( array[ i ] > secondLargest ) {
      if ( i != indexLargestElement )
	secondLargest = array[ i ];
    }
  }
  return secondLargest;
}

/************************************************************************************************************
 * DoNeighborList -- deals with the potential rebuilding of the neighbor list on the basis of the accumulated
 *                   displacement of each site. Adapted from [Rapaport pp.54-56].
 ************************************************************************************************************/
void DoNeighborList( int part )
{
  int n;
  FILE *fPtr;
  real largestDisp, secondLargestDisp;
  int largestDispIndex;

  if ( ( part == 1 ) && rebuildNebrListNow ) {
    /*     fPtr = fopen( "neighborList.log", "a" ); */
    /*     fprintf( fPtr, "step %d -- building neighbor list...\n", stepCount ); */
    /*     fclose( fPtr ); */
    BuildNeighborList();
    rebuildNebrListNow = 0; /* reset Boolean variable */
    DO_SITE dispBetweenNebrListRefreshes[ n ] = 0; /* reset displacements' array */
  } else { /* part 2 */
    DO_SITE {
      dispBetweenNebrListRefreshes[ n ] += ( VLen( siteVelocity[ n ] ) * deltaT );
    }
    largestDisp = FindLargest( &largestDispIndex, dispBetweenNebrListRefreshes, n );
       printf("\nstep %3d: largestDisp = %f\n",stepCount, largestDisp);
    secondLargestDisp = FindSecondLargest( largestDispIndex, dispBetweenNebrListRefreshes, n );
       printf("step %3d: secondLargestDisp = %f\n\n",stepCount,secondLargestDisp);
    if ( ( largestDisp + secondLargestDisp ) > rNebrShell )
      rebuildNebrListNow = 1;
  }
}

